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\  ABSTRACT 

\The  rate  of  convergence  of  line  search  algorithms  based  on  general 

interpolating  functions  is  derived,  and  is  shown  to  be  independent  of  the 

particular  interpolating  function  used.  This  result  holds  for  the  root 

finding  problem  f(x)  =  0  as  well.  We  show  how  inverse  interpolation  can 

be  used  in  conjunction  with  the  line  search  problem,  and  derive  its  rate  of 

convergence.  Our  analysis  suggests  that  one-point  line  search  algorithms 

(in  particular  Newton's  method)  are  inefficient  in  a  sense.  Two-point 

algorithms  using  rational  interpolating  functions  are  recommended. i 
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1 .  INTRODUCTION 


An  essential  part  of  multidimensional  minimization  algorithms  is 
a  line  search,  i.e.,  a  one-dimensional  scheme  for  the  solution  of  the 
equation 

(1)  f'(x)  =  0. 


Most  of  the  line  search  algorithms  in  common  use  are  based  on 
polynomial  interpolation  of  f.  At  iteration  i,  a  polynomial  P  (x)  (the 

XI  f  o 

so-called  hyperosculatory  interpolation  polynomial)  which  coincides  with 
f  and  its  derivatives  up  to  order  s-1 ,  at  each  of  the  n+1  interpolation 
points  x^,  x^  1 ,  xi_n,  is  constructed.  The  new  interpolation  point 
xi+.j,  is  the  solution  of 


(2>  Pn,S(xi+1>=0- 


In  fact,  to  facilitate  the  solution  of  (2),  a  low  degree  polynomial 
is  fitted,  i.e.,  r  =  s(n+1)  is  small;  quadratic  and  cubic  fit  being  most 
commonly  used. 

In  recent  years,  the  possibility  of  using  nonpolynomial  interpolation 
functions  received  some  attention.  One  important  situation  arises  in  line 
searches  associated  with  n-dimensional  constrained  problems,  solved  by 
barrier  function  methods.  A  fit  by  a  polynomial  cannot  capture  the 
singular  behavior  of  the  barrier  objective  function  at  the  boundary  of  the 
feasible  region.  Wright  [20]  dealt  with  the  case  of  the  logarithmic  barrier 
function.  She  suggests  using  the  interpolating  functions 


(3)  ax  +  b  +  r  log(x-c) 

o 

(4)  ax  +  bx  +  c  +  r  log  (x-d). 
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Bj^rstad  and  Nocedal  [3]  analyze  the  rate  of  convergence  of  an 
algorithm  cased  on  the  interpolating  function 


ax  +  bx  +  c 
(dx  +  I)2 


This  function  is  the  one-dimensional  restriction  of  the  "conic''  model 
function  suggested  by  Davidon  [5],  who  lists  some  important  advantages  of 
the  conic  model  over  the  quadratic  one. 


Independently,  we  suggested  [1]  another  rational  interpolating 


function 


ax  +  bx  +  c 
dx  -  1 


which  we  analyze  in  section  3. 

Nonpolynomial  interpolation  was  suggested  much  earlier  for  the  root 
finding  problem 


(7)  f(x)  =  0. 


Ostrowski  [13,  p.  82]  used  in  this  conjunction  the  rational  function  > 

which  Jarratt  and  Nudds  lb]  and  Jarratt  [9]  generalized  to 


x  -  a 
Q(x) 


where  Q(x)  is  a  polynomial.  Ben-Tal  and  Ben-Israel  [2]  describe  nonpolynomial 

interpolations  by  certain  types  of  generalized  convex  functions. 

We  formally  define  the  -  interpolation  algorithm  as  follows. 

n  1 5 

Let  n>0,  s>1  be  fixed  integers  and  let  T  be  a  family  of  s-1  times 
differentiable  functions  T:R+R,  depending  on  r  =  s(n+1)  parameters.  At 


iteration  i,  the  points  x^,  x^,  x^  are  given,  and  a  function 
T  e  ¥  is  chosen  so  as  to  satisfy  the  interpolation  equations 

(9)  T{k)(xi_j)  =  f(k)(xi-j)  *  =  °>  k  =  °»  8~1’ 

A  new  interpolation  point  is  computed  from 

(10)  T'(x1+1)=0  , 

and  the  oldest  point  x,  is  deleted. 

l-n 

The  practicality  of  using  a  particular  class  ¥  depends  to  a  great 
extent  on  the  degree  of  difficulty  of  solving  the  (generally  nonlinear)  system 
of  equations  (9)  and  equation  (10).  In  the  case  of  the  logarithmic  functions 
(3)  and  (4),  equation  (10)  is  easy  to  solve,  but  (9)  is  an  ill-conditioned 
nonlinear  system  of  equations.  Wright  [20]  uses  table  look-ups,  and  applies 
Newton’s  method  after  operating  some  transformations  on  these  equations,  in 
order  to  solve  them. 

For  the  conic  model  studied  by  BjizSrstad  and  Noiedal  [3],  equations  (9) 
and  (10)  can  be  reduced  to  quadratic  equations,  while  for  the  rational  function 
(6)  discussed  in  section  3,  equations  (9)  are  reduced  to  a  linear  system,  and 
(10)  is  very  easy  to  solve. 

Note  that  in  the  polynomial  case,  (9)  is  a  linear  system.  However, 

(10)  is  difficult  to  solve  unless  T  is  a  low  degree  polynomial.  It  will  be 
shown  in  section  4,  that  this  difficulty  can  be  circumvented  by  employing 
inverse  interpolation. 

Note  also,  that  for  the  function  (8),  the  interpolation  equations  can 
be  reduced  to  a  linear  system,  while  the  solution  of  T(xi+1)  =  0  is  simply 
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Investigate  the  rate  of  convergence  of  these 
Here  we  say  that  the  rate  of  convergence  of  a 
to  a  is  p,  if  there  exists  a  positive  number  C, 

C  , 

(see  l 19,  pp.  1-131).  Ortega  and  Rheinboldt  [ll]  refer  to  the  rate  p  defined 
above  as  the  G-order  of  the  sequence  {x^}.  When  it  exists,  it  coincides  with 

their  so-called  Q-  and  R-  orders  (see  [ll,  section  9]). 

Rate  of  convergence  analysis  is  supplied  by  BjtJrstad  and  Nocedal  [3] 
for  the  conic  function  with  s  =  2,  n  =  1.  The  derivation,  which  uses  a 
symbol  manipulation  computer  program,  is  quite  elaborate.  Moreover,  the 
analysis  does  not  carry  over  naturally  to  the  study  of  the  convergence 
properties  of  an  algorithm  using  the  same  interpolation  function,  but  with 
different  data  say  s  =  1,  n  =  3. 

Wright  [20]  gives  no  rate  of  convergence  analysis  for  the  algorithms 
using  the  logarithmic  interpolating  functions  (3)  and  (4). 

An  outline  of  the  paper  is  as  follows.  In  section  2  we  prove  rate  of 
convergence  theorems  for  general  Tn  g  -  interpolation  methods .  We  show  that  the 
rate  of  convergence  is  given  by  the  unique  positive  root  of  the  irdicial 
equation 

(11)  tn+1  -  (s-l)tn  -  s  ^  t^  =  0  . 

J=0 

Since  this  equation  depends  on  n  and  s  only,  the  rate  is  independent  of  the  class  V. 


In  this  paper  we 
minimization  algorithms, 
sequence  {x^}  converging 
such  that 

h+i  ■  al 
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In  section  3  we  analyze  the  specific  family  of  interpolating  functions 
2 

(6)  ax  tbxtc  . 
dx  -  1 

Inverse  interpolation  for  minimization  algorithms  is  introduced  in 
section  4.  We  show  that  the  rate  of  convergence  in  this  case  is  again  given 
as  the  positive  root  of  (11). 

Numerical  examples  illustrating  the  convergence  theorems  are  given 
in  section  5. 

In  section  6  we  discuss  the  implications  of  the  rate  of  convergence 
analysis  to  the  design  of  algorithms. 
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2.  RATES  OF  CONVERGENCE  OF  NONPOLYNOMIAL  ALGORITHMS 

Traub  [19]  studied  the  rate  of  convergence  of  algorithms  that  use 
polynomials  to  interpolate  f,  or  its  inverse  function  for  the  root  finding 
problem  (7).  The  natural  modification  of  these  results  for  the  minimization 
problem  are  discussed  by  Tamir  [l7,  18]  for  the  direct  polynomial  case,  i.e., 
when  the  interpolation  requirements  are  given  by  (9),  T  being  a  polynomial 
of  degree  <  r  =  s(rn-l). 

The  key  result  for  this  analysis  is  the  product  form  formula  of  the 
error  incurred  in  hyperosculatory  polynomial  interpolation  (e.g.  [6,  p.  67]). 
Ostrovski  [l3,  p.  12]  generalized  this  formula  to  the  case  where  the 
interpolating  function  is  not  necessarily  a  polynomial.  However,  no  use  of 
this  generalized  formula  has  been  made  to  extend  the  analysis  of  Traub  and 
Tamir  to  the  nonpolynomial  case.  Using  this  formula,  we  will  obtain  a 
difference  equation  which  differs  from  the  one  obtained  by  Tamir  in  its 
right  hand  side  only.  This  implies  that  in  the  nonpolynomial  case  too, 
the  rate  is  given  by  the  positive  root  of  the  indicial  equation  (11). 

Tamir  [l7,  18]  gives  two  separate  proofs  for  the  cases  s  =  1,  s  >  1.  We 
will  give  a  unified  proof,  and  settle  his  conjectures  in  [17]. 

Stronger  results  than  ours  can  evidently  be  obtained  by  relaxing 
some  of  our  assumptions  (compare  for  example  Brent  [4]  ).  We  have  preferred, 
however,  to  keep  the  presentation  unobscured  by  these  technicalities.  For 
the  same  reason,  we  have  not  stated  explicitly  the  interval  of  (local) 
convergence.  This  is  done  in  great  detail  in  [17]  and  repeated  in  [18]. 

We  will  denote  by  a  a  solution  of  (1),  and  by  J  the  interval 
J  =  (x:  |x  -  a |  <  L}  for  some  positive  L.  The  error  x^  -  a  will  be  denoted 
by  ek,  and  the  open  interval  detexmined  by  {a^,  ...,  affl}  will  be  denoted  by 
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<an,  a^. 


The  following  assumption  will  be  used  repeatedly. 


Assumption  1  r  =  s(n+l)  >3;  f  and  T  have  continuous  derivatives  of 
order  r  +  1  in  J  for  all  Te¥;  f"(a)  i-  0  ;  xi  e  J  and  x^  /  x^  for 
j  ^  k,  i,  j,  k  =  0,  1,  2,  . . . ,  n  ;  e^  i-  0  for  all  k. 

Note  that  if  e^  =  0  for  some  i,  x^  is  a  solution  of  (1),  and  the 
algorithm  is  terminated. 

In  order  that  the  sequence  {x^  defined  by  the  algorithm  be  well 
defined,  the  interpolation  equations  (9) >  as  well  as  equation  (10)  for 
xi+1>  must  have  solutions.  If  T  is  the  class  PQ  g  of  polynomials  of  degree 
less  than  r  =  s(n+l),  equations  (9)  have  a  solution  if  and  only  if 
x^  3*  x^  for  k  4  l.  To  quote  Davis  [6,  p.  27],  the  hope  that  an  interpolation 
problem  can  always  be  solved  providing  the  number  of  parameters  equals 
the  number  of  conditions,  is  naive.  T  can  be  replaced  by  P_  in  iterations 
at  which  (9)  has  no  solution,  but  in  practice  this  case  is  rather  unlikely. 

We  will  assume  henceforth  that  (9)  has  a  solution  for  all  i. 

As  for  equation  (10),  we  will  prove  that  under  Assumption  1,  it  has 
a  solution  for  all  i,  if  L  is  small  enough.  We  need  the  following  difference 
relation  to  prove  this  and  other  results. 

Theorem  1  Under  Assumption  l,  if  T"  ^  0  on  J,  then  the  errors  -  x.,  -  a, 
Induced  by  the  T  -  interpolation  algor  thm,  satisfy  the  recursion  equation: 

I1)S 


n  .  n  n 

w  v'  s-1  _  s  yj-  tmm  s 

Vl  =  1  S  6i-k  j=o  ^  +  1  5=0  9l-J 

j?k 
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where 


M,  = 


Mi(*i<yi))  (-1>r-1  • s  ,  "i  K+i>)  (-x>r  , 


T"  (e(xi+i) ) 


ij.  = 

i 


T"r"e(xi+i)) 


M1(x)  = 


f(r)(x)  -  T(r)(xl 
r! 


N, (x)  =  f(r+1)(x)  -  T(r+1)(x)  , 

1  (r+1) ! 


^(t),  n^(t)  e<t,  x^  x^  y  . xi_n  >  and 
©(Xj^^)  e  ^  >  • 

Proof:  The  error  in  the  interpolation  (9)  is  given  by  (see  [13,  p.  12]  ) 

(13)  T(t)  =  f(t)  +  nQ  (t  -  x^)8  . 

Differentiating  (13)  we  have 

(1*0  f'(t)  =  T'(t)  +  M1(q(t))w(t)  +  N1(ni(t))  W(t)  , 

n  s 

where  W(t)  =  IT  (t  -  x  .)  and  M, ,  N, ,  5, (t)  and  n^(t)  are  defined  above 
j=0  J  ill  i 

(for  proof  see  [l,  section  41  where  we  generalize  Ralston’s  result  [l4,  15]  on  the 
differentiation  of  the  error  term,  to  the  hyperosculatory  case).  Substituting 
t  =  a  in  (14)  and  using 


T'(a)  =  T'(a)  -  T’ (x1+1)  =  -ei+1  T"(©(xi+1))  we  obtain  (12) 


'1+1' 


a 
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Under  Assumption  1,  fn(°0  ^  0.  Since  f ' (a)  =  0,  f 1  must  change  its  sign  at 
a.  It  follows  by  substituting  t  =  a  -  L  and  t  =  a  +  L  in  (14),  that  T*  also 
has  opposite  signs  at  these  two  points  (for  a  detailed  proof  see  Appendix  A 
in  Il8]),  if  L  is  small  enough.  We  summarize  this  result  in 

Theorem  2  Under  Assumption  1,  if  L  is  small  enough,  there  exists  xi+1  e  J 
satisfying  equation  (10). 


Using  Theorem  1  in  [7,  chapter  6,  section  5l>  it  follows  immediately 
from  the  difference  equation  (12),  that  if  the  initial  errors  e0,  ...»  Sq  are 
small  enough  (i.e.,  L  is  small  enough),  the  sequence  tends  to  zero,  establishing 
the  following  local  convergence  result. 


Theorem  3  Under  Assumption  1,  if  L  is  small  enough,  the  sequence  (x^  converges 
to  the  solution  “  of  (1).  □ 

Also  note  that  if  L  is  small  enough,  and  if  s  >  1,  we  have  by  (12) 
lei+ll  <  1  ei'l  >  ^Plying  xi+1  ?  For  s  =  1  however,  we  have  to  assume 

xi+i  t  xi  (cf.  [l6]). 

We  now  replace  (12)  by  a  more  useful  difference  equation. 


Theorem  4  Under  the  assumptions  of  Theorem  3,  and  if  t  0,  then 
(15) 


_  .  s-1  n  s 

ei+l  “  Ai+lei  i-3 


Proof. 


where  Ai+1 — >>H. 

By  assumption,  the  sequences  M.^,  Ni  are  bounded. 


If  s  >  2,  (12)  implies 


(16) 


10 


(i.e.  superlinear  convergence).  If  s  =  1,  we  must  have  n  >2  since  we  assumed 
r  =  s(n+l)  >  3*  For  n  =  2,  (12)  is  the  basic  difference  relation  governing 
the  behavior  of  the  Quadratic  Fit  algorithm,  which  is  known  to  converge 
superlinearly  (see  Theorem  3.4.1  in  Brent  [4]).  It  is  evident  from  (12)  that 
the  rate  for  n  >  2  is  not  less  than  the  rate  for  n  =  2.  Therefore,  (1 6)  holds 
for  all  s  >  1,  n  >  0  if  r  =  Cn+l)s  >  3.  Rewriting  (12)  in  the  form 

n  e,  N, 


U7)  eul  -  Ul} 


1  + 


k=l 


+  M  ei 
i-k  i 


we  see  by  (16)  that  (15)  holds  with 

n 


Aw  *  Mi 


IL  e,  N 

— —  +  -r. —  e. 
e,  M.  i 


k=l 


i-k 


Remark  In  (17,  Appendix  C),  Tamir  conjectures  that  his  apriory  assumption 
[17,  Assumption  2]  on  the  superlinear  convergence  of  the  sequence  {e^}  is 
redundant.  Our  proof  shows  that  this  assumption  is  indeed  redundant. 

We  now  state  our  main  result. 

Theorem  5  Under  the  assumptions  of  Theorem  4,  the  sequence  {x^}  generated 

by  the  T  -  interpolation  algorithm  converges  to  the  solution  a  of  (1),  with 
n,s 

rata  of  convergence  p  which  is  the  unique  positive  root  of  the  equation 

tn+1  -  (s-l)tn  -  s£V  =  0  . 

j=0 

Proof  Convergence  of  (x^  to  a  is  proved  in  Theorem  3.  The  assertion  -Voout 
the  rate  of  convergence  is  a  direct  consequence  of  the  difference  equation  (15), 
as  proved  by  Tamir  [l7,  18]. 

a 
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Remark  The  assumption  that  the  sequence  tL  has  a  nonzero  li  lit,  enables  us 
to  use  the  analysis  of  Tamir  in  [17,  18  J  (or  Traub  in  [  19 3  for  the  root  finding 
problem).  It  will  hold  if  the  mapping  from  ,  •••,  x.,  to  the  parameters  of 
T  defined  by  (9)  is  continuous  so  that  the  limit  exists,  and  if  this  limit  is 
different  from  zero.  In  this  case  the  C-,  Q-  and  R-  rates  of  convergence  are 
exactly  p,  where  p  is  given  in  Theorem  5-  If  the  limit  of  M  exists  and  is 
zero,  or  if  this  limit  does  not  exist,  but  the  sequence  is  bounded,  equation 
(12)  implies  that  the  Q-  and  R-  rates  of  convergence  are  still  at  least  p. 

Corollary  The  rate  of  convergence  of  the  sequence  generated  by  the  interpolation 
algorithm  does  not  depend  od  the  class  of  interpolating  functions  T. 

Remark  It  is  evident  from  our  analysis  that  the  above  corollary  holds  for 
the  root  finding  problem,  as  well  as  for  the  case  when  the  number  of  pieces  of 
information  used  at  the  interpolation  point  x^  depends  on  j  (e.g.  the  False 
Position  Method). 

It  follows  from  Theorem  5,  that  the  rates  of  convergence  of  the  interpola¬ 
tion  algorithms  using  the  conic  interpolating  function  (5)  is  p  =  1.46  for  s  =  .1, 
n  =  3  (**  interpolation  points  with  no  derivatives);  p  =  2  for  s  =  2,  n  =1 
(f  and  f'  used  at  two  points)  and  p  =  3  for  s  =  4,  n  =  0  (f,  f',  f",  f,M  used 

at  one  interpolation  point).  Rates  of  convergence  of  algorithms  using  the 

interpolating  functions  mentioned  in  the  introduction,  can  be  computed  likewise. 

The  behavior  of  the  rate  p  as  a  function  ofn,  for  fixed  s,  is  summarized 
in  Theorem  6. 

Theorem  6  For  fixed  s,  p  is  an  increasing  function  of  n.  For  n  =  0,  p  =  s  -  1, 

while  for  n=i,p  =  s.  As  n  tends  to  infinity,  p  tends  to  + 


Pro°f  For  n  =  0,  n  -  1,  the  rate  is  obtained  by  solving  the  indicial  equations 
t  -  ( s-l)  =  0  and  t  -  (s-l)t  -  s  =  0  respectively.  The  remaining  assertions 
are  proved  in  Tamir  [17,  18].  O 


A  few  numerical  values  for  p,  are  listed  in  Table  2.1. 


TABLE  2.1 


Since  ^  + 


4)  +1  is  close  to  s  even  for  small  values  of  s 


(see  Table  1),  Theorem  6  implies  that  algorithms  using  more  than  two  interpolation 
points  (n  >  1)  are  inefficient.  However,  two  points  algorithms  are  substantially 
faster  than  one  point  (or  memoryless)  algorithms.  Instead  of  making  the  last 
statement  precise  by  defining  a  measure  of  efficiency  (chosen  carefully  to  suit 
the  authors'  purpose),  we  will  note  that  the  transition  from  n  =  0  to  n  =  1 
involves  storage  (but  no  computation)  of  s  extra  pieces  of  data.  In  addition  to 
this,  the  system  of  equations  (9)  will  involve  2s  instead  of  s  unknowns. 

However,  this  system  is  linear  in  the  polynomial  and  rational  cases  (which  are 
the  most  important  ones)  and  need  to  be  solved  once  only  for  the  class  T.  The 
main  difficulty  is  the  solution  of  equation  (10).  This,  in  the  case  of  s  =  3, 
n  =  1  (Newton's  Method  with  memory)  with  polynomial  interpolation,  is  a  polyno¬ 
mial  equation  of  degree  4.  Solution  of  this  equation  can  be  avoided  by  using 

inverse  interpolation,  to  be  discussed  in  section  4.  On  the  other  hand,  for  line 

(k) 

search  algorithms,  computation  of  f  (t)  involves  in  fact  computation  of  the 
derivatives  of  a  function  on  Rn  (i.e.,  gradient  vectors  and  Hessian  matrices,) 
making  the  extra  effort  worthwhile. 
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3.  A  CLASS  OF  RATIONAL  INTERPOLATING  FUNCTIONS 

In  this  section  we  briefly  discuss  the  four  parameter  rational  interpolating 
function 


(18) 


R(x)  = 


ax  +  bx  +  c 
dx-1 


Writing  (18)  in  the  form 

(19)  (dx-l)R(x)  =  ax^  +  bx  +  c  , 


differentiating  (19)  implicitly  and  then  using  the  interpolation  equations  (9), 
leads  to  a  linear  system  of  equations  for  the  coefficients  a,b,c,d.  For  ex¬ 
ample,  with  data  s  =  4,  n  =  0,  the  equations  are 

(dx^-l)f(x^)  =  ax^  +  bx^  +  c 
(dx^-l)f ' (x^)  +  df(x^)  =  2ax^  +  b 
(dXj-Df’^)  +  2df'(xi)  =  2a 
(dxi-l)f"'  (x±)  +  3df'(xi)  =  0  . 


Note  that  if  d = 0,  R(x)  has  no  singularity.  Therefore,  it  may  be  expected 
that  R(x)  will  provide  a  good  fit  to  functions  with  regular  or  singular  behavior. 

We  now  turn  our  attention  to  the  solution  of  (10)  for  xi+^  .  If  d =  0,  R(x) 
is  a  quadratic  and  (10)  yields  x^+^  =  - .  For  d  ^  0,  it  is  convenient  to  re¬ 
write  R(x)  in  the  form 


(20) 


R(x)  =  ax  +  8  + 


15. 


where  a  =  p  =  j  y  =  j  +  -^  +  -~}  g  =  i 

d  d  dZ  d  d2  d3  d 


Differentiating  (20)  we  have 


(21) 


R'(x)  =  a - ; 

(x-&)‘ 


(22) 


R"(x)  = 


iL 


(x-8  )* 


From  (22)  we  see  that  R"(x)  has  exactly  one  change  of  sign  at  x=“6.  The 
point  xi+^  will  be  a  minimum  of  f  if 


(23) 


R'(xi+l)  =  0 


(24) 


R”(xi+i>  >  0  . 


From  (21) -(24)  we  have 


(25) 


xi+1  =  6  ±  •JvJa 


assuming 


(26) 


ay  >  0  . 


The  two  solutions  In  (25)  correspond  to  the  minimum  point  of  the  convex 
branch  of  R,  and  the  maximum  point  of  the  concave  branch  of  R.  Multiplying 
(22)  by  (x-a)4,  we  see  that  in  order  for  (24)  to  hold,  we  must  have  Y(x^+^-&)  >  0, 
which  combined  with  (25)  yields 


16. 


xi+1  =  5  +  e  'ly/a  ,  e  =  sign  y  . 

Condition  (26)  will  hold  near  the  solution  under  the  assumptions  of  Theorem  2. 

Remarks.  Rational  interpolations  are  particularly  useful  in  cases  where  f, 
or  its  derivatives,  have  rapid  changes,  even  when  f  has  no  singularities  (see 
section  5) . 

Use  of  rational  functions  other  than  (5)  and  (6)  suggests  itself,  especially 
when  higher  degree  interpolation  is  needed,  possibly  combined  with  inverse  interp¬ 
olation  (see  section  4). 


* 


t 
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4.  INVERSE  INTERPOLATION  FOR  LINE  SEARCH 

Inverse  interpolation  methods  for  the  root  finding  problem  f(x)=0  are 
well  known.  Assuming  that  f’  is  nonzero  and  f^r)(x)  is  continuous  on  an 
interval  J  mapped  by  f  onto  K,  then  f  has  an  inverse  F,  and  F^r^ 
is  continuous  on  K.  If  T  is  a  hyperosculatory  interpolating  function  satis¬ 
fying 


(27) 


(k}  k 

T  (yi-j)  =  F  (yi-j^  j  =  0,...,n;  k=0 . s  -  1, 


then 


(28) 


•„(r) 


F  (t)  =  T(t)  + 


F  '  (9  (t)  )  -T  '  '  (6  (t)  )  n 


(r), 


r  I 


11  » 
j=0  J 


with  9i(t)  e  (t^  *  yj.i . yi-n)’  In  the  inverse  interpolation  algorithm 

for  the  root  finding  problem,  we  approximate  a  =  F(0)  by  xi+1  =  T(0). 

The  derivatives  of  the  inverse  function  F  can  be  expressed  in  terms  of  the 
derivatives  of  f.  Indeed,  letting 

”k  F  *  °k  ~  *  *  k  ”  1,2 . 

we  have  (see  [12]) 


(29) 


n-k  -1  (2n-k  -2) !  -(2n-k.-l)  k0  k 

Sixhrr*!  1  -2— 

£■  j  n 


where  the  summation  is  taken  over  all  ^  ,  k2  , . . .,  kR  satisfying 
n  n 

S  k.  =  n  -  1  ,  S  ik  =  2n  -  2  ,  k. 
i=l  1  i=l  1  1 


>  0  . 
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Let  T  be  a  polynomial  Q  of  degree  <  r.  By  the  above  and  since 

s 

y.  .  =  f(x.  .)  and  F(y,  ,)=x.  .  ,  Q  can  be  expressed  in  terms  of  the 
i-J  r— J  i"j  n,  s 

(k) 

data  x  .  and  f  <x,  .)•  If  T  is  not  a  polynomial,  we  first  construct 

i'j 

the  interpolating  polynomial  Q  satisfying  (27),  and  proceed  to  solve  the 

n  y  s 


system 


T<'k)(yl-j)  ■  j  =  k  =  0 . 3-1  , 


yi-j  "  f(Xi-j)  ' 


Traub  [19]  shows  that  the  rate  of  convergence  of  the  polynomial  inverse 

interpolation  algorithm  is  given  by  the  positive  root  of  the  (root  finding) 

indicial  equation  t  -  s  £  tJ  =  0,  exactly  as  in  the  case  of  direct  polynomial 

j=0 

interpolation.  Similar  to  our  derivation  in  section  2,  it  can  be  shown  that  the 

rate  of  convergence  is  independent  of  the  interpolating  class  of  functions. 

Inverse  interpolation  has  not  been  applied  so  far  to  the  solution  of  line 

search  problems.  We  will  define  the  ^  -inverse  interpolation  algorithm,  and 

H)  S 

prove  that  under  the  appropriate  assumptions,  its  rate  of  convergence  is  given 
by  the  positive  solution  of  the  indicial  equation  (11). 

A  difficulty  in  applying  inverse  interpolation  to  the  line  search  problem 
is  that  one  cannot  assume  that  f  has  an  inverse  near  an  extremum  point  a, 
since  necessarily  f'(a)  =  0.  Denoting,  however,  g=f'  we  can  write  equation 
(1)  as  g(a)  =  0.  Assuming  that  a  is  a  simple  zero  of  g,  g  has  an  inverse 
G  defined  on  a  neighborhood  of  g(a).  Since  the  solution  a  of  (1)  satisfies 
g(a)  =0,  it  is  given  by 


(30) 


a  =  G(0)  . 


19. 


The  assumption  on  the  differentiability  of  g  implies  that  G  is  differentiable. 
Hence  G  is  continuous  and  has  a  primitive  function  F  (i.e.,  an  indefinite  in¬ 
tegral  of  G),  satisfying 


(31) 


F’(t)  =  G(t)  . 


Equation  (31)  determines  F  up  to  an  additive  constant.  By  (30)  a  is  given  in 
terms  of  any  solution  F  of  (31)  by 

(32)  a  =  F' (0)  . 


Now  let  F  be  any  solution  of  (31),  and  let  T  be  a  hyperosculatory  interpolating 
function  satisfying 


(33) 


T^k\yi_j)  =  j  *  0 . .  k*0,  ...,a  -1  , 

n-i  ■  • 


The  inverse  interpolation  process  for  the  solution  of  (1)  consists  of  approximating 
a  in  (32)  by 

(34)  xi+1=T'(0). 


Evidently,  x^+^  as  defined  by  (34)  is  independent  of  the  particular  integration 

constant  associated  with  F.  Let  Q  be  the  interpolation  polynomial  of  degree 

n,  s 

<  r  satisfying  (33).  We  will  later  express  (34)  in  this  case  in  terms  of  the  data. 
If  T  is  not  a  polynomial,  we  can  express  equations  (33)  in  terms  of  the  data  by 
first  constructing  Q  (i.e.,  replace  T  by  Q  in  (33))  and  then  interpolate 

H) S  Hj 8 

Q  by  T  (i.e.,  replace  F  by  Q  in  (33)). 
n*  s  n,  s 

In  order  to  write  (34)  explicitly  in  the  polynomial  case,  we  can  proceed  to 
construct  P  (x),  the  direct  interpolating  polynomial  determined  by  (9),  differ¬ 


entiate  it  to  obtain 


20. 


n,  8 


r-2 

(x)  =  Z  a,  (x-x  ) 
k=0  K  1 


from  which  we  obtain  directly  by  [12]  the  inverse  interpolation  formula  for 
line  search 


(35) 


i+1 


=  x. 


r-2 

+  z 

k=l 


Pk(“V 


where  is  given  in  terms  of  the  a^'s  in  (29). 

As  an  example,  we  now  construct  the  algorithm  (35)  for  the  case  n~l,  s  =  3. 

This  is  the  algorithm  which  uses  Newton’s  method  data  f,  f'  and  f",  at  the 
interpolation  points  and  x^  ^  .  The  rate  of  convergence  of  this  algorithm  is  3. 

Replacing  x  ^ and  x^  Dy  x3 ,  x2  and  x1  .  respectively,  and  de¬ 
noting  f(k)(Xl)  by  f^k\  we  have 

(36)  P(x)  =  f2  +  (x-x2)f2  +  ^(x-x2)f'2  +  a(x-x2)3 

3  3  2 

+  b(x-x2)  (x-x^  +  c(x-x2)  (x-xL)  , 


where  the  coefficients  a,b,c  are  determined  by  (9)  as: 

^  (f1-f2)-(x1-x2)f2  -  %(x1-x2)2f” 

(xrx2) 

(f {-f 2 ) - (xx -x2) q -3a  (xx-x2) 2 
(x1-x2)3 

(f’^-f'p-6a(x1-x2)-6b(x1-x2)2 
2(x1-x2)3 


i* 

t 


Rearranging  (36)  and  differentiating,  we  obtain 


22. 


Remark .  If  P(x)  is  a  quadratic  (which  is  the  case  for  the  classical  Newton, 
False  Position  and  Quadratic  Fit  methods),  P'(x)  is  a  linear  function,  with 
linear  inverse,  so  that  in  this  case  the  direct  and  inverse  interpolation  formulas 
coincide. 

The  inverse  interpolation  formula  for  s  =  4,  n  =  0  (with  rate  3),  will  differ 
by  the  above  argument  from  the  direct  interpolation  formula  for  this  case.  It  is 
given  by 


(42) 


i+1 


1 

2  3 

(f'P 


Ml 

i 


Note  that  omitting  the  term  with  the  third  order  derivative  in  (42)  yields 
Newton's  method.  Note  also  that  the  direct  interpolation  formula  in  this  case  is 
given  by  the  solution  of  the  quadratic  equation 


P0,3<’1i+l)  *  °>  ”here 

PQ  3(x)  =  fi  +  (x-xi)f|  +  %(x-xi)2f’^  +  !‘(x-xi)3fj"  . 


We  now  turn  to  the  analysis  of  rate  of  convergence  of  this  class  of  algorithms, 
starting  with  the  derivation  of  a  basic  difference  equation. 

Theorem  7.  Let  f "  ^  0  and  let  f^r+i"\  T^r+^  ^e  continuous  on  an  interval  J. 

Let  the  derivative  G  of  F  be  the  inverse  of  g=f',  and  let 

x.  ,,  ,  x.  , . . .,  x.  e  J,  where  x...  =T'(0)  and  T  satisfies  (33).  Then 
i+1  i  i-n  i+1 


(43) 


'i+1 


n 

s 


k=0 


n 


n 

J=0 

j# 


+  L. 


s 

i-j 


J 


where 


The  following  theorem  characterizes  the  behavior  of  the  inverse  interpolatory 
process  for  the  line  search  problem  (1). 

Theorem  8.  Let  f  and  T  have  continuous  derivatives  of  order  r  + 1,  in  the 
interval  J  =  {x:  jx-aj  <L}.  Let  f”(a)^0.  If  L  is  small  enough,  if 
Xq,...,x  e  J  and  the  sequence  (x^)  is  constructed  by  the  inverse  interpolation 
algorithm  for  line  search  (i.e.  x^+^  =  T'(0),  vdiere  T  satisfies  (27)),  then 

x^+^  e  J.  Furthermore,  if  the  algorithm  does  not  terminate,  and  the  sequence 
defined  in  Theorem  7  satisfies  K  4-  0,  then  x^  >  a  with  rate  of  convergence 

p  which  is  the  unique  positive  root  of 

,  i  n-1 

(47)  t  -  (s-l)tn-s  E  t3  =  0  . 

j“0 

Proof.  The  proof  is  identical  with  the  proof  of  Theorem  5. 

Indeed,  both  Theorems  8  and  5  are  consequences  of  the  same  basic  difference 
equation,  namely  (43)  or  (12),  where  the  coefficients  and  L^  are  bounded 

and  lim  ,  lim  are  nonzero. 

Equation  (47)  is  identical  of  course  with  equation  (11)  which  is  the  indicial 
equation  of  the  derived  difference  equation  (15). 
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5 .  NUMERICAL  EXAMPLES 

The  purpose  of  this  section  is  to  illustrate  that  the  theoretical  rate  of 
convergence  predicted  by  the  preceding  theorems,  is  well  reflected  in  the  actual 
behavior  of  the  various  direct  and  inverse  algorithms.  Ten  algorithms  (without 
safeguards)  are  applied  to  minimizing  two  functions: 

f  (x)  =  g*  x^  -  x^  +  2x 

and 

f  (X)  =  X  +  . 

e  -1 

The  first  function,  although  nonsingular,  behaves  very  much  like  a  singular 

one  in  the  interval  [0,2],  due  to  rapid  changes  of  f  and  its  derivatives  in 

the  interval.  This  in  particular  caused  the  cubic  fit  method  to  diverge. 

The  second  function  is  highly  singular  at  x  =  1.  For  this  function,  three 

of  the  methods  based  on  polynomial  interpolation  diverged.  In  contrast,  all  four 

methods  based  on  rational  interpolation  worked  well. 

The  results  are  summarized  in  Tables  5.1,  and  5.2.  The  Rational  and  Conic 

2  2 

functions  referred  to  in  these  tables  are  ~ -  and  - respectively. 

ax-t  (dx+1) 

Initial  values  used  in  Table  5.1  are  (2),  {2,2.1},  {2, 2. 1,2. 2},  and  {2, 2.1, 2.3, 2.3} 
according  to  the  number  of  interpolation  points.  Initial  values  used  in  Table  5.2 
are  {1.75},  {1.7, 1.8},  {1.7,1.75,1.8},  and  {1.7,1.73,1.77,1.8}. 

As  can  be  seen  from  these  tables,  there  is  an  excellent  agreement  between  the 
predicted  rate  of  convergence,  and  the  actual  behavior  of  the  error  sequence  {x^-x  }. 


I 


B 


TABLE  5.1:  Solution  of  f'(x)  =  0  , 


f  (x)  =  g-  x  -  x  +  2x 


26. 


Rational 


Data:  f  at  4  points 


Rate:  1.4 


Polynomial 

(Newton) 


Data:  f,  f 1 ,  f"  at 
1  point 

Rate :  2 


Algorithm 

Iterations 

Polynomial 

No. 

_ 

X 

f'(x) 

* 

X  -  X 

(Quadratic  Fit) 

i _ 

01 

2.1 00000000^0  1 

2.96U01000EJ 

9 . 7925 7 38 8 7'"-' 

"1 

1  1 

:l  .673653.14  U-  0  1 

6 . 72854695  1 1'  0 

5,529105296'-' 

"1 

Data:  f  at  3  points 

-  91 

l  ♦  60/1 2903. l':'0  1 

4 . 972863922l”:0 

4,843864193®' 

"1 

31 

1 » 5 0233685 5 ";0  1 

2. 882030 163K0 

3, 8 15942  4  4  OK' 

"1 

41 

1 .39553441)8^0  1 

1.346951822'^) 

2,44791.8768® 

"1 

Rate:  1.3 

.  :;i 

.1  .3  1  036/301'“  0  1 

7.684444291  """'1 

1 .976246894'-" 

"1 

61 

1  .25724  1530®0  1 

3.99216/669®  "1 

l. 364989 189® 

"1 

71 

1 . 2095  1  0432'"-0  1 

1  , 997530860':-""  1 

8 . 876782049®" 

* . 

81 

I.I/A099919K0  1 

1 .00566836 7® ""1 

5 . 5  35  73  0723'- 

91 

L .  151783312*0  1 

4. 7.1  86 3 54 4 OK ""2 

3.104070026® 

-9 

A- 

101 

.1 .135973  L55r-!o  1 

2. 03 4 304 IV UK "2 

1. 523054 197® 

t  1.  1 

1  .  l76917/08'-  0  1 

7.6.1  3347509'" ‘  3 

6 . 1 75096 25 OK 

"3 

1  9  1 

1  ,  1226003  16'- 0  1 

2 . 210373360''  '  3 

.1.  ,865705073'- 

_3. 

Rational 


Data:  f,  f  at  2  points 


Rate:  2 


oi  ioooooooo'-o 

II  I. .  49/9/8645'"0 
21  I  .331855539'  0 
.51  1,  J  94055957'0 

41  1.14  1 25 I 46  1 1  0 

hi  I  .  I  2279.{VM6"  0 
A  I  I  .  I  2'>  7606V  I  •'  0 

n  i .  120/424151  <> 


01 

2.300000000**0 

1 

5 . 049343000® 1  1 

1 . 1 792 57389® 0 

l  1 

1.548461 1  70""  0 

1 

3, 70909.1  269*0  1 

4.2771. 85466®'“ 

21 

1. .  463037249"  0 

1 

2,281  68484 1'  0  1 

3.422946381®'“ 

.31 

1 . 352784382*0 

1 

1.04 03 895 4 OK 0  1 

2 , 3  2 0 4 1 7 7 1 0 ®  ” 

4  1 

1  .25611 4652'- 0 

1 

3 , 93661. 3400 '"""I.  1 

1 « 3537204  1 1'5""" 

51 

1..  1  8  708666 6 "O 

1 

1  . 2°  76 096 6 3 '••"1  1 

4.634405477®'"" 

61 

1  . 1534.179761-0 

1 

5.030661 751 K- 2  1 

3,26753642  L'""' 

71 

l  ,  .1  32 6 87573 1 0 

1 

1.5  5  06  3  1.81  IK  “2  | 

1 . 194496.1 1  7''""" 

81 

l .  12359  4309"- 0 

1 

3. 4090304 78®"  3  1 

2,8516977  1. 3'~" 

91 

l .  12 1.09  7655®  0 

1 

4.  146773605'"- "”4  1 

3 , 550432350'' "• 

.1  01 

1 .  120760307'-  0 

1 

2 . 0601 97  l38r-"5  1 

1 . 76 95 78 524 '•""" 

t  1  1 

1 .  12074281  7' "0 

1 

2 . 394566790K" 7  1 

2,056193661'"“ 

121 

t ,  1  207426  1  2''  0 

1 

3.6  1. 6562856'" 10  1 

2, 3.66502078'""“ 

01 

2 , 000000000^0 

1 

2. 20  00  OOOOOi- 1  1 

8 . 7  9 2 5  7 3 8 8  7 "  1 

1.  1 

l . 6 76 4 70508® 0 

1 

6 , 811 13522 8'"  0  1 

5 « 557279769'-""! 

21 

1 .  4 4 5092 362"' 0 

1 

2.0371 18820*0  1 

3 . 2 4349 7 5 1 1 '“ ""  .1. 

31 

1 . 289992758'  0 

1 

5 , 7 9 9  6 094  6 2 1  1 

1  . 69250 1468'"'"":l 

4  1 

1 ,  L  950086 98 '  -0 

1 

1 .528615287"'""  l  1 

7  ♦  426608709':""2 

51 

1.14450  1.369*0 

1 

3. 407897060®  "2  | 

2 . 375875747'*“2 

6 1 

1  .  124595007'  0 

1 

4.649413295":  "3  1 

3 . 852395469'“  ""3 

71 

1 .  120075207''  0 

1 

1 .54633845 7" "4  1 

1  .326760935''  "4 

!.)'  1 

1  .  120742778'  0 

1 

l  ♦  945599880'-'  "7  1 

.1.  .670493971'-'  7 

91 

L .  120  74261 J  i-O 

1 

3.09  4629587""'  131 

“9 , 3777, 68484"  "  1.1 

2 .  96  I  I  01 000 ".I 
..1 . 0  I  0902098";0 
8,691  73  l  487'  "  .1 
t  .5,'39950i2'-"  l 
2 .  HA  4  I  4505  I  "  ‘  2 
2 ,  4353/96  |  Ji-  3 
3. 036 53 7-5  7vr 
4,41 1  'obh/yi  v 


9, 7925 7300 /K"  l 
3.  7 7 2 .5 6 0.5 34* “l 
2.1  LL 129277'--  1 
7  .411  334555®  *:> 

2.0500049/8'.:  "2 

2 ,  Ob 20. .-5 4  22  41 "  3 
2.  4 07905933®  '5 
3,695909565'  "9 


o 


TABLE  5.1  Continued 
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Algorithm 


Iterations 


Inverse  Polynomial 


No. 


f’(x) 


X  -  X 


01  2.1 00000000^0 


Data: 

f, 

f*  at  2  points 

1  1 

1 . 620084808^0 

1 

5.4869441 38^0 

1 

6.280848878"=" 

1 

21 

1 . 1 9050873 5"=0 

1 

1 . 395258096*'“ 

1  1 

1.905087348*" 

1. 

.  31 

7, 014376:142^'“  1. 

1 

6.937507796'”' 

1  1 

" 2.98 5 6238 5  8 "= " 

1 

Rate: 

2 

41 

9.706749253'='“:!. 

1 

1.185756774'-"“ 

1  I 

"7 . 932507474"=" 

2 

51 

9.808130  l.40'=~.1. 

1 

2.  169400795'“'' 

7  1 

'“1 . 918698516"=" 

9 

61 

9 .  yyi  j  1 2J  09'"“  j 

1 

8.9431 16557'"" 

4  1 

"8. 88789051 2'=“ 

4 

71 

9.999902581'“  “1 

1 

1  ,741909068'"" 

4  1 

“l. 741887828 

6 

81 

1,0000000001=0 

1 

6, 889045:1  07'"  "' 

121 

“'6,889045541"" 

'12 

Conic 

01 

2 . 10000000  O' "0 

1 

2,961.  10  1.000'-' 1 

1 

9 , 792573887"='" 

'1 

1  1 

1 .45835  1778f:0 

1 

2. 2.1  409091  6*0  1 

3.376086170"=" 

"1 

.  21 

1 .296 13036  I  P  0 

1 

6. 10 1.350832'"“ 

1  1 

1  ♦  753877493"='" 

1 

Data: 

f, 

f'  at  2  points 

,51 

1. 1.72934975*0 

1 

9 , 776009320'- " 

7  1 

5 , 2194-34386'“" 

"2 

4  1 

1 .13  161 22 70 tO 

1 

1.397730063'"'“ 

2  1 

1 . 086966669"='" 

“2 

Rate: 

51 

.!  .  .121  263007' 0 

1 

6. 00756442  1"”" 

4  1 

5.203952432'='" 

"4 

2 

6  1 

l  .  120744056" 0 

1 

1 . 48  1  494620'  - " 

6  1 

1  , 444434108":" 

"6 

7  1 

1.12074261  l'  0 

1 

1  , 1  26  799  608*“' 

1  1.  1 

“8 . 436336612"=' 

"11 

_  • 

01 

2.  100000000'“0 

1 

2 ,941 1.0  1. 000*1 

1 

• 

9 . 792573887"=" 

"1 

inverse  ruiynuuu.aj, 

.  1  1 

1 .497496006' 0 

1 

2.803094043^0  1 

3.767533944'"" 

"1 

Data: 

f. 

f',f" 

71 

1.216743480''') 

1 

2.236092670*'" 

1  1 

9 . 550086866"=' 

-9 

Am 

at 

2  points 

31 

1 . 178643239'--) 

1 

9. 8 94  7 098  7  3' r 

3  1 

7 . 900627235ft" 

"3 

Rate: 

O 

4  1 

1  .  1  70/44.1  39'=  0 

I 

L .  778781  894*"" 

4  1 

.1. .  528009852"=“ 

6 

J 

51 

.1  .  L 207 4261  1 ' - 0 

1 

8,6/36  1.7380'"" 

I9| 

"'9 . 404354075"=' 

"u 

2.96U01000E1 


1.100000000E0 


Inverse  Polynomial 

Data:  £,f’,  at 

1  point 

Rate:  3 


0  1 

7.000000000"".) 

1 

2 . 700000000"=  l 

1 

8. 792573887'=  “1 

1 1 

1 ,55/945/56"'  0 

1 

3 . 894  70267  4 *=0 

1 

4 , 372031448"=  "1 

21 

1.299.1  50/28""  0 

1 

4,3/4399616"-"  .1 

1 

1 .  /84031 1 62*''"  1. 

3  1 

.1. .  1714.1  424  4  0 

1 

8. 95/4580  to*""' 2 

1 

5  087163538'"  "2 

4  1 

L. 1264840 19*0 

1 

7.053276/68* "3 

1 

5. 74340787  l'=  “3 

51 

.1 .  17.0/48788'- 0 

1 

3,04 //80  175"=  "5 

1 

2 , 61 7642078"=  “5 

41 

t .  120/4761  l'=0 

1 

3. 55765  1647"='  "'  1 

21 

""9. 0991508  L6'“  “11 

Rational 

Data:  at 

1  point 

Rate:  3 


01  2 . 000 0000 00*=0 
'.II  l  ,4/2542767"'=0 
7!  1 . 277734985* 0 

I  .  134799/10"0 
41  1  .  I  2O8H430 /"  0 

lil  1. .  120 /'1 2 6.1  L'-O 


7. 200 0 000 00i- 1 
2.4J  854 7443'=  0 
2.4/8962995'  "I 
I  ,858405152'  "  "7 
1 . 45  I  599  944'='  •  A 

2 .  3  0  0  /  2  2 5  29"““  10 


8 . 792573887*'"  1 
3 .  ?  j  1. 7  9  <?  <f>  “7 15  B  •  =  “  1 
1 .0,1  9923737ft"  1 
1 ,405709866'=- "2 
1 .4 1699402/'“  ""4 

1 .036082452"="  10 
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TABLE  5.2:  Solution  of  f'(x)  =  0  ,  f(x)  =  x  +  l/(ex"*- 


1) 


Algorithm 


Iterations 


Polynomial 
(Quadratic  Fit) 

Data:  f  at  3  points 


No. 


f'(x) 


Rate:  1.3 


01  1 . 800000000^0 
’Ll  1.8981*80/11-0 
21  1.933113075E0 

■31  1.953688313*0 

41  1.961024208*0 

51  L  ♦  96226251 .1  *0 

6  1  1.962  4 1.6421 '-0 

7  I  1 . 96 24 23527'- 0 

81  1.962423649*0 


"4.81 7670946* "1  I 
"1  . 59562  4554*"' 1  I 
'6. 86741. 3860*'  2  I 
1.98021  446  A-"  2  I 
"2.1261 16528*"3  I 
3. 60409801 4'-“  “4  I 
"L. 61642 4827''"  '5  I 
'2.762306271*'  7  I 
2 . 808237059'“  ”9  I 


"1, 62423650 1*""1 
"6 . 4 2 5 5 5 7 8 7 5 rr ” 2 
"2 .93 1057467* "2 
”8.735336800E~3' 
'1. '399442557'"  "3 
1. 6.1. 1395:1  09*'~  4 
'7. 228  789794'"  “6 
“1  .2353764571"-' "7 
"1 , 25592651  le~9 


Rational 

01 

1.800000000*0 

1 

'"  4 , 8 1 7 67094  6 *  “ '  1 

1  1 

1 .9623434681-0 

1 

"1.793161898* "4 

Data:  f  at  4  points 

21 

31 

1 . 962403344k  0 

1 . 962423642*0 

1 

1 

'4 ,093437  L  4  4  *"’5 
"  1 . 80 6 0 9693 2 E" 8 

Rate:  1.4 

41 

L . 962423650*0 

1 

't.  46548 7887E-1 

"  1 .624  2  3 6  5 0 1.  e '"  .1 
"8, 018257350 K ”5 
"1.8305882901^  "'5 
'8. 077 11 0926* ”9 

'  6  ♦  5 5 3  8  6  3  3 5 5 *  "13 


Polynomial 

(Newton) 

Data:  at  one 

point 
Rate:  2 


01 

1. .  7 50000000'-' 0 

"6. 96736890 l*"!  | 

'"2. 12423650  L* 

"1 

1  1 

1 ♦ 897153528*0 

"1.62 3 305337 E"1  1 

'"6,527012224'" 

"2 

21 

1.95591 2208*0 

'"1 .47 096 5849 f""' 2  1 

"'6 « 51 1  392281'" 

”3 

31 

1 .962357443^0 

'.1  ♦  400!i949U8':r  "4  1 

" 6 .62073590 4* 

"5 

4  1 

1 .962423643':  0 

"1.53  4  L  58280*' "8  1 

"6.860964333* 

"9 

51 

.1 . 962423650'-0 

"1.6479B'7302'':'"'16I 

'"  "7. 3  7  257  4  773  r:. 

"17 

Rational 

Data:  f,f'  at  2  points 
Rate :  2 


01 

1.800000000*0  1 

1  '"4.8!  76 7 09 4 6* "l  1 

"1  ♦  62423650 1*' 

"  .1 

1  1 

1 .9623432031:0  I 

1  "1 .7543 132  74E~'  4  1 

"7 . 844698287*' 

"5 

21 

1  .962 4 2363 9* 0  1 

1  ~ 2,4417 8 9574*"'  8  1 

'"1  .092001476'"' 

"8 

31 

1 » 962423650*0  1 

"1 .86482  7737*"'  1  71 

"8,23993651.1  *' 

~lfl 

Conic 

Data:  f,  f'  at  2  points 
Rate:  2 


01  1*800000000*0 
II  1 . 969742383'^) 
21  1 .962364402*0 

31  1.962423651*0 


"4.81  76)70946'  “l  I 
L  .61  79591 93* '"2  I 
'1.3249514 26E“4  I 
1.8627.L2947'"~9  I 


M.  62423650  L  i"  "1 
7. 31873279 l  fe— 3 
"5.92481341  I H-5 
8. 3 3030 555 6 1C 


Inverse  Polynomial 
Data: 

at  one  point 
Rate:  3 


01  1 ,750000000''0 

11  1.940533164*0 

21  L  .962395  4 93'  -0 
31  t . 962423650^0 


"6. 96 736890 IE  "1  | 
5 .06780 9 7 6 5 f;; '”  2  I 
'6. 296  301  78 1'-""  5  I 
"1.364  3  5 5  6  7  7 * ""  1  3  I 


"2. 124236501 '“  "1 

'2.189040593' . 2 

"2  ♦  8  L  5703434'“  "5 
6 . 1 0 1  5  8  6> 2 5  0 1  ■  “  1 4 


Rational 

Data: 

at  one  point 
Rate:  3 


01 
l  I 


1  .  75 00000 OOf-  0 
1 . 96234  I 099*0 
1  . 9624  2 3 65 0&:  0 


6.967368901'-'"  l  I 
"t  .046146248*"" 4  I 
"1.397579968':  '"1 4  I 


"2 . 12 4236501 1"  "  L 
0,255 1502 1 4*""S 
"  6. 2502086  04  K"  If: 


6.  CONCLUDING  REMARKS 


Our  analysis  (and  limited  numerical  experience)  suggest  that  rational  interp¬ 
olating  functions  should  be  preferred  over  polynomials;  the  rate  of  convergence  is 
the  same,  and  for  both  rational  functions  and  polynomials  the  system  (9)  is  linear. 
However,  rational  functions  can  better  cope  with  singularities  and  perform  equally 
well  for  regular  functions. 

The  analysis  also  points  to  the  inefficiency  of  interpolation  algorithms  based 
on  more  than  two  interpolation  points  (or  more  than  three  points  if  function  values 
only  are  used).  Two-point  algorithms  are  significantly  faster  than  one-point  algor¬ 
ithms,  the  latter  are  therefore  useful  only  if  computation  of  the  derivatives  of  f 
are  relatively  very  cheap. 

Use  of  inverse  interpolation  is  recommended  if  equation  (10)  is  difficult  to 
solve.  Note  that  even  in  the  Cubic  Fit  case  where  the  interpolating  function  is  a 
cubic,  solution  of  equation  (10)  involves  computation  of  square  roots  (see  [10,  p. 

142]),  in  itself  a  relatively  costly  operation  on  the  computer  (cf.  equations  (41)). 

Finally,  note  that  any  modifications  made  in  the  algorithms  in  order  to  ensure 
convergence  (see  e.g.  [10  section  7.3]),  may  severely  affect  the  rate  of  convergence, 
since  the  basic  difference  equations  may  be  fundamentally  changed  by  such  modifications. 

Assume,  for  example,  that  we  modify  the  Quadratic  Fit  algorithm  so  that  one  of 
the  points  x^+^  ,  ^  ^  (not  necessarily  x^  ^ )  is  dropped,  in  such  a  manner 

that  the  remaining  points  bracket  the  solution  a.  Then  we  may  choose 
e^  <  0  <  <  e^  and  small  enough  L,  such  that  equation  (15)  of  Theorem  4  would 

imply  that,  for  M  >  0,  we  have  e^+^  >  0  for  all  i.  Hence,  in  the  bracketing 
algorithm,  one  of  the  three  interpolation  points  is  fixed  as  x^  ,  and  in  the 
difference  relation  (15)  one  of  the  indexes  should  be  replaced  by  3,  leading  to 
difference  equation  with  an  indicial  equation  different  than  (11). 


30. 


Thus  the  statement  in  Tamir  [17],  that  bracketing  algorithms  do  not  lend 
themselves  to  the  difference  equation  approach,  and  the  conjecture  made  there 
that  the  interpolation  and  the  bracketing  algorithms  have  the  same  rates  of 
convergence,  are  both  false. 

A  bracketing  procedure  that  aims  at  maintaining  the  rate  of  convergence  of 
the  underlying  interpolation,  should  coincide  with  it  near  the  solution. 
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